************************************************************************************************
************************************************************************************************
*
* Title:  Decay or Resilience: The Long-Term Social Consequences of 
*         Conflict-Related Sexual Violence in Sierra Leone
*
* Author: Carlo Koos
*
* please make sure you have the following add-on programs installed:
*        - coefplot
*        - graph scheme "plotplain"
************************************************************************************************
************************************************************************************************

************************************************************************************************
************************************************************************************************
*
* Preparations
* 
************************************************************************************************
************************************************************************************************
* save the do-file and the dataset in the same folder
* load dataset <change to directory where data and do file are located>
clear
cd "<ENTER DIRECTORY>"
use decay_resilience_sierra_leone.dta, clear

* Run Coarsened exact matching (CEM) procedure to balance the data and reduce model 
* and specification dependence
imb HQ_RUF hhsize fatheredu poverty urban muslim , treatment(raped)
cem HQ_RUF(0 1 3 7) hhsize(3 7 11) fatheredu(0 1 3 4) poverty urban muslim , treatment(raped)

* surveyset with new_weight
* new_weight combines both the household data calculated weight of an observation 
*   and the cem-weight from the CEM procedure
* wta_hh = sampling weight for the primary sampling unit from the household survey
* cem_weight = weight of one obs after cem procedure
g new_weight=wta_hh*cem_weight

* survey set with new_weight
svyset, clear
svyset slihseacode [pweight=new_weight], strata(stratum)

************************************************************************************************
************************************************************************************************
*
* Set covariates for models
* 
************************************************************************************************
************************************************************************************************
global X1 "raped                                                                 urban hhsize poverty muslim female age age2"
global X2 "raped                                hhmkilled limbs displaced HQ_RUF urban hhsize poverty muslim female age age2"
global X3 "raped  humaid                   hhmkilled limbs displaced HQ_RUF urban hhsize poverty muslim female age age2"
global X4 "raped  humaid rapedXhumaid hhmkilled limbs displaced HQ_RUF urban hhsize poverty muslim female age age2"

************************************************************************************************
*
* Table 1: SUMMARY STATISTICS
* 
************************************************************************************************
eststo clear
svy: logit memcommunity $X3 if cem_matched==1 , noconstant  // this model runs only to calc the descriptive stats for the sample
estpost sum memcommunity sqrtsocialcont sqrtdonations $X3 if e(sample) 
esttab using descriptive.rtf, replace cells("mean(fmt(%9.3f)) sd(fmt(%9.3f)) min(fmt(%9.0f)) max(fmt(%9.0f)) count(fmt(%9.0f))") noobs label ///
   title("Table 1: Summary Statistics") onecell compress
    
************************************************************************************************
*
* Table 2: EFFECT OF CRSV ON MEMBERSHIP IN ASSOCIATIONS
* 
* Note that I use STATA’s survey module to account for the first stage (district level) and the 
* second stage, i.e. the primary sampling unit (PSU), meaning that all models use standard 
* errors clustered at the PSU.
*
* Note: I keep the unmatched sample in order to run robustness checks further below. In the 
*       main analysis I only use the matched sample.
*
************************************************************************************************
eststo clear
svy: logit memcommunity $X1  i.district if cem_matched==1 
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store A

eststo clear
svy: logit memcommunity $X2  i.district if cem_matched==1 
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store B

svy: logit memcommunity $X3  i.district if cem_matched==1 
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store C

svy: logit memcommunity $X4  i.district if cem_matched==1 
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store D

*output for main text (without control vars)
esttab A B C D using association.rtf, ///
  replace label compress nogaps drop(_cons) ///
  scalars("N Number of observations" "fstat F-Stat" "fpval p-Value" "model Model") ///
  title(Effect of CRSV on Membership in Associations) order($X4) b(3) se(3) obslast ///
  indicate("Household controls = hhsize* poverty* urban*" "Individual controls = muslim* female* age* age2*" "District dummies = *.district") ///
  star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
  addnotes("\textit{Note:} Outcome variable \textit{association} measures whether a household member is member of a community-based association. Robust standard errors clustered by primary sampling unit in parentheses.")
  
*output for appendix (with control vars) APPENDIX IN TEX FORMAT
esttab A B C D using reg1full.tex, ///
  replace label eqlabels(none) booktabs compress nogaps drop(_cons) ///
  scalars("N Number of observations" "fstat F-Stat" "fpval p-Value" "model Model") ///
  title(Effect of CRSV on Membership in Associations \label{mainfull}) order($X4) b(3) se(3) obslast ///
  indicate("District dummies = *.district") ///
  star(+ 0.10 * 0.05 ** 0.01 *** 0.001) nonotes ///
  addnotes("\textit{Note:} Outcome variable \textit{association} measures whether someone in the household is member of" ///
    "a community-based association." "Robust standard errors clustered by primary sampling unit in parentheses." ///
	"+ \textit{p} < 0.10, * \textit{p} < 0.05, ** \textit{p} < 0.01, *** \textit{p} < 0.001")

************************************************************************************************
*
* Table 2: EFFECT OF CRSV ON PROSOCIAL BEHAVIOR
* 
* Note that I use STATA’s survey module to account for the first stage (district level) and the 
* second stage, i.e. the primary sampling unit (PSU), meaning that all models use standard 
* errors clustered at the PSU.
*
************************************************************************************************
eststo clear
svy: logit memcommunity $X4 i.district if cem_matched==1
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store A
margins, dydx(*) post
est store AM

svy: logit sqrtsocialcont $X4 i.district if cem_matched==1
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store B
margins, dydx(*) post
est store BM

svy: logit sqrtdonations $X4 i.district if cem_matched==1 
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store C
margins, dydx(*) post
est store CM

*output for main text (without control vars)
esttab A B C using prosocial.rtf, ///
  replace label compress nogaps drop(_cons) ///
  scalars("N Number of observations" "fstat F-Stat" "fpval p-Value" "model Model") ///
  title(Effect of CRSV on Membership on Prosocial Behavior) order($X4) b(3) se(3) obslast ///
  indicate("Household controls = hhsize* poverty* urban*" "Individual controls = muslim* female* age* age2*" "District dummies = *.district") ///
  star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
  addnotes("\textit{Note:} Outcome variable \textit{association} measures whether a household member is member of a community-based association. Robust standard errors clustered by primary sampling unit in parentheses.")

*output for appendix (with controls)
esttab A B C using regalloutcomesfull.tex, ///
  replace label eqlabels(none) booktabs compress nogaps noconstant drop(_cons )  ///  
  scalars("N Number of observations" "fstat F-Stat" "fpval p-Value"  "model Model") ///
  title(Effect of CRSV on Prosocial Behavior \label{alloutcomesfull}) order($X4) b(3) se(3) obslast ///
  indicate("District dummies = *.district") ///
  star(+ 0.10 * 0.05 ** 0.01 *** 0.001) nonotes ///
  addnotes("\textit{Note:} Outcome variable \textit{association} measures whether someone in the household is member of" ///
    "a community-based association." "Robust standard errors clustered by primary sampling unit in parentheses." ///
	"+ \textit{p} < 0.10, * \textit{p} < 0.05, ** \textit{p} < 0.01, *** \textit{p} < 0.001")

************************************************************************************************
*
* Figure 3: EFFECT OF CRSV ON COMMUNITY RELIANCE
* Based on marginal effects from previous models AM, BM, CM
*
************************************************************************************************
coefplot AM, bylabel(association) || BM, bylabel(contributions) || CM, bylabel(donations) ||, ///
  scheme(plotplain) level(90) byopts(row(1)) xline(0)  grid(none) drop(*.district rapedXhumaid baseline urban hhsize poverty muslim female age age2) // play(avgmarginal)

// change the plot graphically  
gr_edit .plotregion1.subtitle[1].style.editstyle drawbox(no) editcopy
gr_edit .plotregion1.subtitle[1].style.editstyle size(small) editcopy

graph export marginal_effects_three_models.png, replace width(16000)

************************************************************************************************
*
* Figure 4: EFFECT OF CRSV ON COMMUNITY RELIANCE
* The underlying regression model is also estimated here, but the table is only shown in the appendix.
*
************************************************************************************************
label define raped 0 "CRSV in HH=0" 1 "CRSV in HH=1"
label values raped raped

eststo clear
svy: logit assist_community raped hhmkilled limbs displaced HQ_RUF urban hhsize poverty ///
  muslim female age age2 i.district if cem_matched==1 
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store A

margins, at(raped=(0 1)) atmeans 
marginsplot, recast(scatter) scheme(plotplain) allxlabels level(90) title("") ytitle("Pr(Reliance on Community)") play(margins_plot)
graph save reliance, replace 
graph export reliance.png, replace width(16000) 
  
*regression table output for appendix
esttab A using reliance.tex, ///
  replace label eqlabels(none) booktabs compress nogaps noconstant drop(_cons humaid)  ///  
  scalars("N Number of observations" "fstat F-Stat" "fpval p-Value"  "model Model") ///
  title(Social Acceptance \label{reliance}) order($X3) b(3) se(3) obslast ///
  indicate("Household controls = hhsize* poverty* urban*" "Individual controls = muslim* female* age* age2*" "District dummies = *.district") ///
  note("Robust standard errors clustered by primary sampling unit in parentheses.") star(+ 0.10 * 0.05 ** 0.01 *** 0.001)

************************************************************************************************
*
* Table 4: FACTORS ASSOCIATED WITH CRSV
* 
************************************************************************************************
global selection "hhsize girlratio fatheredu poverty urban muslim HQ_RUF"
eststo clear
eststo: svy: logit raped $selection if cem_matched==1
esttab using selection.rtf, ///
  replace label compress nogaps drop(_cons) ///
  scalars("N Number of observations") ///
  title("Factors associated with CRSV") b(3) se(3) obslast ///
  star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
  addnotes("Robust standard errors clustered by primary sampling in parentheses.")
    
************************************************************************************************
*
* Section "PARTICULARITIES, LIMITATIONS, AND ADDITIONAL TESTS"
* Percentage of CRSV-affected households among (1) ever displaced/never displaced households and 
*     (2) those who returned/did not return to their home communities
* 
************************************************************************************************
tab raped displaced if cem_matched==1, col
tab raped returned if cem_matched==1, col


************************************************************************************************
************************************************************************************************
************************************************************************************************
*
* APPENDIX
*
************************************************************************************************
************************************************************************************************
************************************************************************************************

************************************************************************************************
*
* Figure 1: INTERACTION EFFECT CRSV IN HH*HUMANITARIAN AID (in the appendix)
*
************************************************************************************************
global XXX "i.raped##i.humaid  i.hhmkilled i.limbs i.displaced c.HQ_RUF i.urban c.hhsize c.poverty i.muslim i.female c.age c.age2"
* int effect for association
svy: logit memcommunity $XXX if cem_matched==1 , noconstant  // cluster(slihseacode)
margins raped, at(humaid=(0 1)) post noestimcheck
marginsplot, xdimension(at(humaid)) level(90) recast(scatter) recastci(rarea) scheme(plotplain)  ///
  title("") ytitle("Pr(Association)") xtitle("Support")  ///
  legend(pos(6) rows(1)) plot1opts(lpattern("")) plot2opts(lpattern("--"))
graph save rape_support_membership.gph, replace
* int effect for contribution
svy: logit sqrtsocialcont $XXX if cem_matched==1 , noconstant  // cluster(slihseacode)
margins raped, at(humaid=(0 1)) post noestimcheck
marginsplot, xdimension(at(humaid)) level(90) recast(scatter) recastci(rarea) scheme(plotplain)  ///
  title("") ytitle("Pr(Contribution)") xtitle("Support")  ///
  legend(pos(6) rows(1)) plot1opts(lpattern("")) plot2opts(lpattern("--"))
graph save rape_support_contribution.gph, replace
* int effect for association
svy: logit sqrtdonations $XXX if cem_matched==1 , noconstant  // cluster(slihseacode)
margins raped, at(humaid=(0 1)) post noestimcheck
marginsplot, xdimension(at(humaid)) level(90) recast(scatter) recastci(rarea) scheme(plotplain)  ///
  title("") ytitle("Pr(Donation)") xtitle("Support")  ///
  legend(pos(6) rows(1)) plot1opts(lpattern("")) plot2opts(lpattern("--"))
graph save rape_support_donation.gph, replace

graph combine rape_support_membership.gph rape_support_contribution.gph rape_support_donation.gph, ///
  ycommon r(1) c(3) colfirst iscale(0.8) scheme(plotplain)
graph export rape_support_int.png, replace width(4000) height(3000)

************************************************************************************************
* Tables 1 to 3 in the appendix are the extended tables of the main text including the 
*               coefficients for household and individual control variables
*
************************************************************************************************


************************************************************************************************
* Table 4: Main results for Association with unmatched sample
* (Appendix)
*
************************************************************************************************

svyset, clear
svyset slihseacode [pweight=wta_hh], strata(stratum)

eststo clear
svy: logit memcommunity $X1  i.district  //, cluster(slihseacode)
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store AU
*
eststo clear
svy: logit memcommunity $X2  i.district  //, cluster(slihseacode)
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store BU
*
svy: logit memcommunity $X3  i.district  //, cluster(slihseacode)
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store CU
*
svy: logit memcommunity $X4  i.district  //, cluster(slihseacode)
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store DU

*output for appendix (with control vars)
esttab AU BU CU DU using reg1fullunmatched.tex, ///
  replace label eqlabels(none) booktabs compress nogaps drop(_cons) ///
  scalars("N Number of observations" "fstat F-Stat" "fpval p-Value" "model Model") ///
  title(Effect of CRSV on Membership in Associations (unmatched) \label{reg1outcomesunmatched}) order($X4) b(3) se(3) obslast ///
  indicate("District dummies = *.district") ///
  note("Robust standard errors clustered by primary sampling unit in parentheses.") star(+ 0.10 * 0.05 ** 0.01 *** 0.001)

************************************************************************************************
*
* Table 5: Main results for Prosocial Behavior with unmatched sample
* (Appendix)
*
************************************************************************************************
eststo clear

svy: logit memcommunity $X4 i.district  //, cluster(slihseacode)
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store A

svy: logit sqrtsocialcont $X4 i.district   //, cluster(slihseacode
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store B

svy: logit sqrtdonations $X4 i.district   //, cluster(slihseacode)
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store C

*output for appendix (with controls)
esttab A B C using regalloutcomesfullunmatched.tex, ///
  replace label eqlabels(none) booktabs compress nogaps noconstant drop(_cons )  ///
  scalars("N Number of observations" "fstat F-Stat" "fpval p-Value"  "model Model") ///
  title(Effect of CRSV on Prosocial Behavior (unmatched) \label{alloutcomesfullunmatched}) order($X4) b(3) se(3) obslast ///
  indicate("District dummies = *.district") ///
  note("Robust standard errors clustered by primary sampling unit in parentheses.") star(+ 0.10 * 0.05 ** 0.01 *** 0.001)  
  
************************************************************************************************
*
* ALTERNATIVE OUTCOMES:
* Table 6: Membership in professional associations
* 
************************************************************************************************
label var memprof "Professional association"

eststo clear
svy: logit memprof $X4  i.district if cem_matched==1  //, cluster(slihseacode)
estat gof
estadd scalar fstat = r(F)
estadd scalar fpval = r(p)
estadd local model "Logit", replace
est store memprof

*output for appendix (with control vars)
esttab memprof using outcomealt.tex, ///
  replace label eqlabels(none) booktabs compress nogaps drop(_cons) ///
  scalars("N Number of observations" "fstat F-Stat" "fpval p-Value" "model Model") ///
  title("Effects of CRSV on Membership in Professional Associations" \label{outcomealt}) order($X4) b(3) se(3) obslast ///
  indicate("District dummies = *.district") ///
  note("Robust standard errors clustered by primary sampling unit in parentheses.") star(+ 0.10 * 0.05 ** 0.01 *** 0.001)

************************************************************************************************
*
* Table 7: 2SLS IV Estimation
* Instrumental Variable Results (Appendix)
*
************************************************************************************************
* diaclose: if diamond mine is less than 10km away from survey location
g diaclose=0
replace diaclose =1 if dist_diamond<10 
* Below: This is when a mine is 10k distant max + medium RUF presence is in the area: 
* medium presence proxies the strat to control more area without having the full 
* relation with the local pop
g diaruf=0
g med_RUF=0
replace med_RUF=0
replace med_RUF=1 if HQ_RUF >4 & HQ_RUF<8
replace diaruf=0
replace diaruf=dist_diamond*med_RUF
label var dist_diamond "Dist. from mines"
label var diaruf "RUF and mines"
*gen "inequality" variable to block potential backdoor path at the IV estimation
recode s9q07 (1=5) (2=4) (3=3) (4=2) (5=1), gen(inequality1) // measure for perceived inequality
label var inequality1 "Inequality"

eststo clear
eststo: ivreg2 memcommunity (raped=dist_diamond diaruf) hhmkilled limbs displaced hhsize poverty inequality1 ///
  urban muslim female age age2 [pweight=wta_hh] if cem_matched==1, r first savefirst savefprefix(first_) 
eststo: ivreg2 sqrtsocialcont (raped=dist_diamond diaruf) hhmkilled limbs displaced hhsize poverty inequality1 ///
  urban muslim female age age2 [pweight=wta_hh] if cem_matched==1, r first savefirst savefprefix(second_) 
eststo: ivreg2 sqrtdonations (raped=dist_diamond diaruf) hhmkilled limbs displaced hhsize poverty inequality1 ///
  urban muslim female age age2 [pweight=wta_hh] if cem_matched==1, r first savefirst savefprefix(third_) 

esttab first_* est1 est2 est3 using instrument0full.tex, se replace label booktabs compress nogaps ///
  b(3) se(3) order(dist_diamond diaruf raped) ///
  stat(N widstat idstat idp, labels("N" "F-Stat weak ident." "LM test underident." "Underident. p-value")) ///
  star(+ 0.10 * 0.05 ** 0.01 *** 0.001) drop(_cons) ///
  title(2SLS IV Estimation\label{ivest})  note("Robust standard errors in parentheses.") ///
  mgroups("First stage" "Second stage", pattern(1 1 0 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) ///
  mtitles("DV: CRSV in HH" "DV: Association" "DV: Contribution" "DV: Donations")

*END


